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Abstract 

Experience with core-collapse supernova simulations shows that accurate accounting of total 
particle number and 4-momentum can be a challenge for computational radiative transfer. This 
accurate accounting would be facilitated by the use of particle number and 4-momentum transport 
equations that allow transparent conversion between volume and surface integrals in both configu- 
ration and momentum space. Such conservative formulations of general relativistic kinetic theory 
in multiple spatial dimensions are presented in this paper, and their relevance to core-collapse 
supernova simulations is described. 
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I. INTRODUCTION 



The state of the art in core-collapse supernova simulations now includes energy- and 
angle-dependent neutrino transport 0, ^, |, ||, |[. However, experience shows that 

simultaneous conservation of both energy and lepton number is difficult numerically 
This challenge motivates us to develop new conservative formulations of relativistic kinetic 
theory that are specifically attuned to the need for accurate accounting of particle number 
and energy in numerical simulations of radiative transfer problems. 

Before describing the conservative formulations of kinetic theory we seek, we explain why 
energy- and angle-dependent neutrino transport is necessary in supernova simulations and 
detail the magnitude of the challenge of energy conservation. 

Sophisticated treatments of neutrino transport are necessary because the ultimate en- 
ergy source of the supernova explosion — the gravitational potential energy of the stellar 
progenitor's core — is eventually converted almost completely into neutrinos. Some of the 
gravitational potential energy is lost to escaping neutrinos during core collapse, but most 
of it is converted into a thermal bath of dense nuclear matter, photons, electron/positron 
pairs, and trapped neutrinos deep inside the nascent neutron star. Neutrinos, having the 
weakest interactions, are the most efficient means of cooling; they diffuse outward on a time 
scale of seconds towards a semi-transparent region near the surface of the newly forming 
neutron star, and eventually escape with about 99% of the released gravitational energy. 
In modeling the conversion of gravitational potential energy into neutrino fluxes, energy- 
and angle-dependent neutrino transport is necessary to accurately follow the transition from 
quasi-isotropic diffusion to forward-peaked free streaming. In this transition region energy is 
transferred from the neutrino radiation to the matter behind a stalled shock wave, and this 
energy transfer may be necessary to propel the shock through the outer layers in an explo- 
sion |jlO|, |l^ . But whether or not such neutrino heating is the proximate cause of explosion. 



the fact that neutrinos dominate the energetics implies that accurate neutrino transport is 
integral to any realistic and comprehensive study of the explosion mechanism. 

The importance of accurate neutrino transport is a lesson learned from experience with 



supernova simulations. Parametrized studies |T2[ highlight the sensitivity of explosions to 
neutrino luminosities and to conditions in the semi-transparent region near the nascent 



neutron star's surface (see also Ref. ||T3[) — precisely the region where neutrino energy and 
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angle dependence must be tracked carefully. Moreover, there remains a nagging quantita- 
tive uncertainty in simulations with multidimensional hydrodynamics: Those with neutrino 
transport that depends on direction in configuration space but is averaged over energy and 
angles exhibit explosions |T^, |T^, 0, while those with neutrino transport that depends on 
energy but is averaged over angles in both configuration and momentum space do not show 
explosions |T^, |TB|. It may be that these differing outcomes are due to the different neutrino 
transport schemes, both of which are ultimately inadequate. 

Moving beyond these general arguments for the necessity of accurate neutrino transport, 
quantitative consideration of the energetics shows how severe the requirements are on one as- 
pect of accuracy: energy conservation. As mentioned above, virtually all of the gravitational 
potential energy (~ 10^^ erg) released during collapse is eventually converted into intense 
neutrino fluxes lasting several seconds. However, supernova explosion energies (the kinetic 
energy of the ejecta) are observed to be only ~ 10^^ erg. Now, because it is difficult to argue 
with any rigor about the physical impact of any energy lost or gained due to numerical 
error, the total energy should be conserved to a precision corresponding to the phenomena 
of interest in the problem. Hence a simulation's result for explosion energy accurate to, say, 
10% would require total energy conservation to an accuracy of about 0.1%. Allowing for 
systematic error accrual, the total energy would have to be conserved at a level of 0.1%/A^ 
per time step, where N ~ 10^ is the total number of time steps in the simulation. 

Conservative formulations of kinetic theory would help to meet the numerical challenge 
of particle number and energy conservation in core-collapse supernova simulations. To give 
an idea of the kind of formulation of kinetic theory that we seek, we first review familiar 
descriptions of the dynamics of a fiuid medium. The dynamics can be described in two 
different ways, which might be called elemental and conservative. 

The elemental formulation expresses the evolution of the fluid in terms of equations of 
motion for the velocity and some independent set of quantities (e.g. temperature, densities 
of various species comprising the fluid, etc.) measured by an observer moving along with 
the fluid ("comoving observer"). For example, consider a spacetime with metric components 
{OtMu} and metric determinant g = det{gf^^), containing a perfect fluid with 4- velocity com- 
ponents {m'^} and comoving frame total energy density p, pressure p, and baryon density n. 
lin the absence of radiative transfer and signiflcant energy input from nuclear reactions, the 
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perfect fluid evolves according to 




(1) 



(2) 



(3) 



(Greek and latin letters are spacetime and space indices respectively.) Supplementary re- 
lations between p, p, and n — referred to as the equation of state — are determined by the 
microphysics of the fluid. The name elemental denotes the fact that by writing down sep- 
arate equations of motion for the velocity and comoving-frame quantities, the kinetic and 
"intrinsic" fluid energies — two "elements" of the system — are analytically separated. 

The conservative approach expresses the evolution of the system in terms of the divergence 
of the stress-energy tensor T^'^. For a perfect fluid, T^^= (p + p)u^u^ + pg'^'^, and Eqs. (|I]) 
and d^) are replaced by 



Volume integrals of the left-hand sides of Eqs. (|) and (|^) — obtained by multiplying by 
the invariant spacetime volume element y/—g d'^x and integrating — are related to surface 
integrals in an obvious manner. Physically, this relates the time rates of change of 4- 
momentum and baryon number in a volume to fluxes through a surface surrounding that 



existence, we comment on what it means for conservation issues. 

There are several important cases where the connection coefficients F'^p^ on the right- 
hand side of Eq. (^) might not vanish. They are present in curved spacetime, where (at 
least in part) they embody gravitational forces. But even in flat spacetime, coordinates 
employed by accelerated observers give rise to connection coefficients. And even without 
spacetime curvature or accelerated reference frames, connection coefficients arise from the 
use of curvilinear coordinates. 

What does the right-hand side of Eq. (^ mean for 4-momentum conservation? Only 
when it vanishes — that is, only for inertial observers in flat spacetime employing rectilinear 




(4) 



while Eq. (H) is replaced by 




(5) 



volume; hence the labeling of Eqs. (^) and (^ as conservative. 

The right-hand side of Eq. (^ deserves note. After discussing the reasons for this term's 
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coordinates — are the components of total 4- momentum constant in time ("conserved"). For 
only in this case do the coordinates reflect the translation invariance of flat spacetime, the 
physical origin of 4-momentum conservation. (Curved spacetime lacks translation invariance, 
so there is no 4-momentum conservation.) Because the presence, for whatever reason, of a 
source term like the right-hand side of Eq. (^ means that the 4-momentum components in 
such a basis are not conserved, it might more properly be called a "balance equation" than 
a "conservation law" . But because the volume integral of the left-hand side easily translates 
into a surface integral, we still call it a "conservative" formulation. 

There are some special cases in which a conserved quantity associated with the time coor- 
dinate t can be found, however. For unaccelerated observers in fiat spacetime, r*p^ = 0, even 
in curvilinear coordinates. This means that energy is conserved, though the 3-momentum 
components in cuvilinear coordinates are not. Another special case is the restriction to 
spherical symmetry in general relativity: Here certain coordinate choices allow the non- 
vanishing r*p^ terms to be absorbed into the left-hand aside, leading to the identification of 
a conserved energy-like quantity (see e.g. Ref. ||T9| ). 

Having used the familiar example of a perfect fluid to discuss what we mean by "ele- 
mental" and "conservative" formulations, we now consider kinetic theory in terms of these 
categories. The evolution of a particle type described by kinetic theory is often expressed 
as an equation of motion for the distribution function /, the ensemble- averaged density of 
a given particle type in phase space. (These concepts will be defined with greater precision 
in subsequent sections; for the present discussion it is sufficient to understand that phase 
space is the combination of configuration space and momentum space, and that multiplying 
/ by the volume of an infinitesimal cell in phase space gives the number of particles having 
positions and momenta within the ranges defined by that cell.) The distribution function 
evolves due to advection through phase space and particle interactions. 

Advection through phase space gives rise to derivatives of / with respect to the com- 
ponents of the position vector x and the momentum vector p. For numerical evolution it 
is convenient to parametrize the distribution function in terms of {x'^}, the components of 



X in a global "coordinate basis" [^. (Here and throughout this paper, quantities defined 
with respect to the coordinate basis have indices without accents.) For momentum we make 
a different choice of basis, because interactions with a fiuid are most easily described — and 
best handled numerically — in terms of momentum components measured by a comoving 



5 



observer. The change from the coordinate basis to an orthonormal basis associated with the 
comoving observer (a "non-coordinate basis" ) has two parts. First there is a transformation 
dx^ = e^^dx^ to an (in general non-comoving) orthonormal basis. (Here and throughout 
this paper, quantities defined with respect to the non-comoving orthonormal basis have in- 
dices accented with a bar.) This is followed by a Lorentz transformation dx^ = A^^dx^ to a 
comoving orthonormal basis. (Here and throughout this paper, quantities defined with re- 
spect to the comoving orthonormal basis have indices accented with a hat.) Hence advection 
through phase space will involve derivatives of / with respect to the coordinate basis position 
components {x^} and the comoving orthonormal basis momentum components {p'^}. 

Turning from advection to particle interactions, we consider the case where the particle 
species are sufficiently dilute that the interactions can be described in terms of a collision 
integral C[/] depending only on the distribution functions / of the individual particle species. 
(This approximation ignores correlations between particles; that is, the number of instances 
of finding particles at the same position is obtained from the product of their distribution 



functions.) In this case, the equation of motion for the distribution function / is [21 



/(AV^^-rW|^)=c[/]. (6) 

This is the Boltzmann equation. 

In terms of the categories described above in connection with fluid evolution, the Boltz- 
mann equation is "elemental" : The most fundamental quantity — the distribution function — 
is the evolved variable, and volume integrals of the equation in both configuration space and 
momentum space are not obviously related to surface integrals. This "non-conservative" 
character is present even in fiat spacetime and rectangular coordinates. The first term of 
Eq. (P) is non-conservative even in fiat spacetime and rectilinear coordinates because the 
boost — which depends on the coordinates {x^} — sits outside the derivative df /dx^. 
The spatial dependence of the boost also gives rise to non- vanishing T'^^^, even in fiat space- 
time and rectilinear coordinates; therefore the second term of Eq. (|^) does not vanish. This 
second term is non- conservative because the factor pP sits outside the derivative df/dp^. 

While the Boltzmann equation is "elemental" or "non-conservative", it is well known 



e.g., see Refs. [20, 22, B3]) the first two momentum moments of the distribution function 



(integrated over a suitable invariant momentum space volume element dP) constitute a 
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particle number current A^^ and particle stress-energy tensor T^^", 

= j fP^dP = j fp^K^[,e^-,dP, (7) 
T^"" = J flfp''dP = j fp^A-^f,e^f,p'A%e%dP, (8) 

which (for electrically neutral particles) obey the balance equations 

1 d 



/—g dx'^ 
1 d , 



-g dx^ 



^^m) = j C[f]dP, (9) 

— = -r%^r^^ + J c[fYdP. (10) 



While this result is often stated in the literature (e.g., see Refs. |2^, |2^), because 
of the non-conservative character of Eq. (^) it is not obvious how its momentum moments 
give rise to Eqs. (|^) and (p^O]). Equation (|^) contains factors A'^^e^^ outside df /dx^] but 
from Eqs. (n)~(P^ we see that these factors have come inside the derivative with respect to 
x^. What happens to the spacetime derivatives of A^f^e^p that are generated in taking these 
factors inside the derivative? Furthermore, according to Eqs. (^ and (§), the quantities A^'^ 
and T^^" involve no momentum derivatives of /. But the second term of Eq. @ contains 
a factor pP{df /dp'^). Because of the momentum factor outside the momentum derivative 
it is not obvious how this term will contribute to Eqs. (P) and ( p!OD when integrated over 
momentum space. 

How, then, is the connection between Eq. (|^) and Eqs. (|^) and (|^) established in detail? 
The reviews of relativistic kinetic theory by Lindquist and Israel do not provide 



detailed proofs. As discussed in Sec. |ITT|, an explicit proof by Ehlers relies on the fact 
that A^'^ and T^^ are momentum-integrated quantities, and no direct insight is gained into 
what happens to the momentum derivatives of / in the integration over momentum space. 

For those interested in computer models of radiative transfer problems, these are not 
idle academic questions; they are issues that must be faced in order to build simulations 
capable of making meaningful scientific statements about the core-collapse supernova explo- 
sion mechanism. Experience with supernova simulations in spherical symmetry shows that 
understanding the detailed connection between the Boltzmann equation and the particle 
number and 4-momentum balance equations has important consequences for how well these 
quantities are conserved in a simulation 24]. While a discretization of the Boltzmann 



equation is a natural numerical method of evolving the neutrino species, naive differencings 



of the various terms in Eq. (^) generally will not be consistent with a straightfoward dif- 
ferencing of Eqs. and (p!0|), leading to unacceptable numerical errors in particle number 
and energy conservation. In Lagrangian (or "comoving" ) coordinates in spherical symmetry, 
Mezzacappa and Bruenn derive a conservative formulation of the Boltzmann equation 
transparently related to particle number balance as expressed in Eq. (^). They also devise 
methods of handling momentum derivatives that are consistent with both number and en- 
ergy conservation Liebendorfer et al. P] went a step further in this spherically symmetric 
case, deducing the connnection between the Mezzacappa and Bruenn "number conservative" 
Boltzmann equation in comoving coordinates and energy conservation as represented by an 



Eulerian (or "lab frame") version of Eq. (10). Using complicated, non-intuitive differencings 
of hydrodynamic and gravitational variables, they construct a numerical implementation of 
the Mezzacappa and Bruenn "number conservative" Boltzmann equation that is stable and 
faithful to its analytic connection to the lab frame version of Eq. ([T0|) to the accuracy nec- 
essary to make solid scientific statements about the neutrino-driven explosion mechanism in 
spherical symmetry ^ 

Because spherically symmetric models with Boltzmann transport fail to reproduce some 
important observable characteristics of core-collapse supernovae (e.g. the launch of an ex- 
plosion |T], 1^, 1^, |, ^ ^), this work must be followed up in multiple spatial dimensions (see 
Refs. [0, 1^ for some early efforts). In this paper we develop — allowing for full relativity and 
multiple spatial dimensions — conservative formulations of kinetic theory, such that volume 
integrals in both configuration and momentum space are trivially related to surface inte- 
grals. These conservative expressions make transparent the connection between Eq. (|^) and 
Eqs. (H) and ([l0|). They can be used to deduce the term- by-term cancellations involved in 
this connection, thereby illuminating the complicated differencings required to achieve the 
cancellations numerically. Our conservative formulations also suggest new primary variables 
of radiation transport: particle number and energy variables describing the contribution of 
each comoving orthonormal basis momentum bin to the coordinate basis number and energy 
densities. It may be that the use of these new radiation variables could provide a simpler 
path to an accurate accounting of particle number and energy in simulations of radiative 
transfer problems. 

The organization of this paper is as follows. Differential forms and exterior calculus are 
natural mathematical tools for handling the volume elements and integrations needed in rel- 
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ativistic kinetic theory. We closely follow (and slightly extend) Ehlers' work p2[ in reviewing 
their use in the description of phase space for particles of definite mass (Sec. y) and the 
derivation of the Boltzmann equation (Sec. |111| ). The centerpieces of this paper — two con- 
servative reformulations of the Boltzmann equation, which provide transparent connections 
to particle number and energy-momentum balance as expressed in Eqs. (|^) and (|I^) — are 
presented in Sec. |rV|. Because differential forms and exterior calculus may be unfamiliar to 
those whose primary interests are radiation transport, in general, or supernova science, in 
particular. Sees. P-|rV| each will contain two subsections: one containing a general deriva- 
tion, and a second that explicitly demonstrates aspects of the derivation in the familiar case 
of the 0{v) limit in Lagrangian coordinates in spherical symmetry. (While we review some 
aspects of exterior calculus in our presentation, these are mostly in endnotes, and are more 
in the character of reminders than a self-contained introduction. For the latter we refer 
the reader to Refs. p2|, ^.) A conclusion (Sec. 0) summarizes our results, comments 



on their connection to moment formalisms, and discusses the utility of these formulations 
for supernova simulations. As an application of our formalism, an appendix contains 0{v) 
equations in fiat spacetime, but in coordinates sufficiently general to represent rectangular, 
spherical, and cylindrical coordinate systems. 



II. PHASE SPACE FOR PARTICLES OF DEFINITE MASS 

In this section we consider phase space for particles of definite mass. We specify particle 
trajectories — the "states" whose average occupation is specified by the distribution function 
/ — and the volume elements needed to derive macroscopic observables from / and operators 
that act upon it. In discussing the general case (Subsec. |11 A| ) we follow Ehlers closely, 
and refer the reader to his work for more detailed discussions and proofs of some of the 
assertions made here. We extend his discussion to the use of momentum components mea- 
sured in a frame comoving with the fiuid that interacts with the particle species treated by 
kinetic theory. The specific case of Lagrangian coordinates in spherical symmetry to 0{v) 
is treated in Subsec. |11 B| . 



9 



A. Phase space for particles of definite mass: the general case 



A study of kinetic theory begins with consideration of the trajectories of individual par- 
ticles. The worldline x{X) of a particle of mass m with 4-momentum p and electric charge e, 
moving through a spacetime with metric components {gfj,^} and electromagnetic field tensor 
components {F^i,}, is determined by 

= eF^^p'' - V^pP^pP. (12) 

In the coordinate basis associated with the position vector components {x'^}, the connection 
coefficients are 

■p/i 

2^ \dxP ^ dx- dx-J- ^ ' 

As mentioned in Sec. |, we seek results expressed in terms of (orthonormal) momentum 
components {p^} measured by an observer comoving with the ffuid with which the particles 
interact. Combining the transformation e'^^ to an orthonormal tetrad and boost A^^ to the 
comoving frame into the composite transformation 

C% = A^^e%, (14) 



Eqs. ( pTD and ([T2|) are replaced by 



^'f.P^ (15) 



= eF^^.p- - r^p^pV- (16) 
The connection coefficients do not transform exactly as tensor components, but as 

^ up ^ p^ p-L up ^ ^ p^ p ■ 

This completes the desired specification of particle trajectories. 

The set of position vectors x and momentum vectors p comprise an 8-dimensional manifold 
M, the one-particle phase space for particles of arbitrary rest masses. The positions x are 
points in spacetime, and the momenta p are points in that portion of the tangent space at 
X characterized by < and p future-directed. The curves (x(A),p(A)) obtained from the 
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particle equations of motion comprise the phase flow in M. The phase flow is generated by 
the Liouville vector 

Here ^ {gfTr} , |^| ^ have been chosen as basis vectors on M. 

Specification of a particle mass m defines a hypersurface of M, a 7-dimensional manifold 
called the one-particle phase space for particles of definite mass m. The mass satisifies 

= -g.uP^P" = -g.uC^i.C'uP^p'. (19) 

Considered as a scalar function m{x,p) on M, the particle mass satisfies 

L[m] = 0, (20) 

which expresses the constancy of m on each phase space trajectory. Hence the 7-dimensional 
manifold is generated by all the phase space trajectories of particles of mass m. Equation 
(pop shows that L is tangent to M^, and Eq. (|19]) is a constraint indicating that it is one of 



the momentum dimensions of M that has been lost in going over to M^. Hence the Liouville 
operator restricted to Mm can be expressed 

Lm = + (eFh/ - r^.,pV)|f , (21) 

where ^ {gfir} ' {^} ) have been chosen as basis vectors in M^. To allow for a change of 
coordinates in momentum space from {p*} to {u*} (to momentum space spherical coordi- 
nates, for example), we rewrite Eq. (|2l|) as 



i". = p'^".^ + (-^-'^Z - (22) 

where ^ {gfjr} ' |^| ) chosen basis vectors on Mm- 

Before discussing volume and "surface" elements on Mm, we give a brief reminder of 



some properties of differential forms and their exterior derivatives (e.g., see Refs. ||2^, ^). 
If {dz"'} are basis 1-forms on some manifold, an m-form ip can be expanded as 

i^=—,^a,...a^dz'''-''-, (23) 

ml 

where dz"-^-"-"^ denotes the wedge product dz"-^ A ... A dz""^. The quantities ipai...am cire the 
components of ip- The components are completely antisymmetric {ipa-L.-am = '^[ai...a^])5 as is 
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the wedge product, so that ijj can also be expressed 

V' = ^|a,..a^|Ci^"^-"'", (24) 

where the vertical bars indicate that summation is taken only over Oi < 02 < • • • < Qm- The 
eocterior derivative takes an m-form into an im + l)-form as follows: 



Two properties of the exterior derivative will be used later in this paper: 



c?((/, + -0) = d(l) + di/j, (26) 
d{(j)Ailj) = d(j)Ailj + {-iy(j)Adi/j, (27) 

where is a p-form. 

Now we turn to a discussion of volume and "surface" elements on M^. The invariant 
volume element in spacetime is the 4-form 



V = ^V^e,,padx'''"^, (28) 



where g is the determinant of the metric, e^i,pa is the completely antisymmetric symbol 
with 60123 = +1, and the wedge product of basis 1-forms is abbreviated by the notation 
dx^^'^" = dx^ A dx'' A dx'^ A dx'^ . Shifting attention to the momentum space at a given space- 
time point, an invariant volume element on a mass shell corresponding to mass m in the 
tangent space is 



= ^/J^^oijkdf'. (29) 

Expressed in terms of basis 1-forms conjugate to orthonormal comoving frame momentum 
components, it becomes 

- ^dp'^', (30) 



where the bold character denotes a 3-vector, and E{p) = ^/\p\'^ + w?. After a final change 
of momentum coordinates from p to u, the momentum space volume element can also be 
expressed as 

' ■ (31) 



E(p) 
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A volume element on Mm is a 7-form, which can be expressed as 



1] An, 



Eip) 



det 



[-] 



123 



(32) 
(33) 



The "surface element" normal to the phase flow in is an important quantity; it is a 
6- form, obtained by contracting [^] Q with the Liouville vector Lm'- 



Lrn. ' ^ 



2!E(p) 



det 



[-] 



du'' 



— e 



(rfw^" Ar/). 



dpi 



dikn 



(34) 



The first term in Eq. ( p^]) can be written in terms of surface elements in spacetime. A 

surface element in spacetime with normal in the x'^ direction is given by the contraction of 

the vector d/dx^ with the spacetime volume element rj: 

d 



dxf^ 



V 



-9 



In terms of these spacetime surface elements, the first term in Eq. 



-9 



3! 



(35) 

can be expressed 

(36) 



This first term of the "surface element" u is all that contributes to an integration over 
momentum space |3^; it corresponds to a cell in the familiar 6-dimensional Newtonian 
phase space. The surface element u also shares an important property with a Newtonian 
phase space volume element: It is invariant under phase flow. This result follows from the 
vanishing exterior derivative of u, 

duj = 0, (37) 



and the use of Stokes' theorem on a small tube of worldlines 22 



B. Phase space for particles of definite mass: Lagrangian coordinates in spherical 
symmetry to 0{v) 

Here we specialize certain results of the previous subsection to massless, electrically 
neutral particles propagating through flat spacetime, employing Lagrangian coordinates in 
spherical symmetry and retaining only terms to 0{v). 
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First we introduce these coordinates. In flat spacetime and Eulerian spherical coordinates, 
the metric components are specified by the line element 

ds^ = -dp + dr^ + de^ + P sin^ 9 d(f)'^. (38) 

We seek Lagrangian coordinates t and m to replace the Eulerian coordinates t and r. The 
enclosed baryon mass m is defined by 

m = Aur'^pdr, (39) 



where p is the baryon mass density. The mass density obeys the 0(f) conservation law 
where 

/A-t\ 

" = Ti 

is the radial velocity measured by an Eulerian observer. We choose the Lagrangian time 
coordinate such that the metric components are specified by a line element of the form 

ds^ = -df + gmm dm^ + dO"^ + sin^ 9 (42) 

The transformations between the Eulerian and Lagrangian coordinates involve eight non- 
trivial quantities, the elements of the Jacobian matrices 

(at at \ 
dr r 1 ' 
at am / 

(at m \ 
at dr / 

Two of these elements, dm/dr and dm/dt, are determined by Eqs. (|39|) and (^Ol). The 
relationship between Eqs. (|38|) and (^2D — specifically, the requirements gu = 1 and gtm = 
— provide two more equations. The identity J!^ir}^{t,m} = {J{t,m}~>{t,r}) ^ provides the 
remaining four equations, and we find that to 0{v) 



' at 


at 


at 


am 


ar 


dr 


. at 


am 


at 


at 


at 


dr 


dm 


am 


at 


dr 



1 

V 



1 

1 -V 



-Anr'^pv Airr'^p 
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(45) 
(46) 



From these relations we see that to 0{v), the hne element in our Lagrangian coordinates is 



-di' + ( -A^ ) dm^ + dO'^ + sin^ 9 (47) 



sin^^ ^ 



whose determinant is 



These metric components can be used to obtain results valid to 0{v) 

Now we consider the equations of motion that determine the particle trajectories. 
The rates of change, along a particle worldline, of the coordinate basis position com- 
ponents {x^} = {t, m, 6, (f)} and comoving orthornomal basis momentum components 
IpAj = are given by Eqs. ([T5|) and (p^): 

dx^ 



dX 
dpf" 
~d\ 



C^f^pr (49) 

-r^p//. (50) 



We now consider the quantities C^jx and T^pp appearing in these equations. As discussed in 
the previous subsection, = A'^^e'^^ is a composite transformation from the coordinate 
basis (denoted by unadorned indices) to an orthonormal comoving basis (denoted by hatted 
indices). In the present case the coordinate basis is already comoving, so that the boost 
A^^ is simply the identity transformation. Under the remaining transformation e^^ to an 
orthonormal basis, the metric of Eq. (^) must become the Lorentz metric. We make an 
obvious choice and take the non-zero transformation elements C^^, to be 

C\ = 1, (51) 
1 

47j-r2p' 

C\ = r, (53) 
= rsinO. (54) 

The elements of the inverse transformation are obvious. Turning to the connection 
coefficients, we employ Eqs. (pTS]), (^Tj), and (^Sj) to find the non-zero coordinate basis 



connection coefficients T^^^p to be 
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T\e = rv, (56) 
= rvsm^e, (57) 

r^«.-r-..^-(^ + ^), (58) 

^2TTr^p ~^ dm ) ' ^^^^ 
T'^ee = -ATfr^p, (60) 

r™^^ = -Airr^psin^e, (61) 

r^, = r%, = r^^ = T%t = ^, (62) 
r me = r gm = r'^m</( = '^'^(pm — ^^^s^' (63) 

= -sin^cos^, (64) 

r^<^ = rV = cot^. (65) 

Now we can obtain the connection coefficients in the orthonormal comoving frame by em- 
ploying Eqs. (O), (EH)"® (^5|)-(|65l). We find that the non-zero F'^,)^ are 

rO„.r.„^-(fi.^). (aa) 

po,, _ yo.. - r^.. - r3.. - Y. (a7) 

22 ~ 33 ~ 20 ~ 30 ~ ^ ' \^ ' ) 

r^22 = r^33 = ~r^2i ^ ^^si = (68) 

^9 ^•^ cot 9 

r'33 = -r'32 = (69) 

Note that unhke the V^^p, the r'^j>p are not symmetric in their lower indices up because the 
comoving orthonormal basis is not a coordinate basis. With C^j^ and F'^p^ now explicitly 
specified, our consideration of the particle equations of motion is complete. 

The next topic to consider from the previous subsection is the Liouville vector. This vec- 
tor generates the phase fiow — the set of trajectories {x{X),p{X)) obtained from the particle 
equations of motion, Eqs. ( ^Q]) and (pOl). These are trajectories through the 8-dimensional 
manifold comprised of position vectors x and momentum vectors p. But a 7-dimensional 
submanifold, the phase space for particles of definite mass, is of more immediate interest: 
The physical fact that the particles have definite mass — as embodied in the "mass shell" con- 
straint of Eq. (p]9| ) — implies that there are only three independent momentum variables. We 
choose these to be {m*} = {e,p,(p}, related to the comoving orthonormal basis momentum 



16 



components {p*} by 



= e — /i^ sinyj. 
In the present case of massless particles, the mass shell condition gives 



(70) 
(71) 
(72) 

(73) 



From these expressions one can compute the Jacobian [dp^/du^) and its inverse {du^/dp^); 
this inverse is 



/ de 

dp^ dp^ dp^ 

d/i dji d/j, 

dp^ dp^ dp^ 

y dip^ d(f dip 
dp^ 



dp^ 



\ 



( ^ 



\ — COS a/1 — /U^ sin y9 \ 



^yjl—^^ cos ip 



V 



— sin tp 



cos tp 



(74) 



We are now ready to write down the Liouville vector on the phase space for particles of 
definite mass, obtained with the help of Eqs. (H), (il)-(il), and (|70D-(|74D: 



9,-2 d 2 

Lm = e— + Arcr pe/i- h e 

at om 



^ /3f ^ 91np 
^ \ r dt 



d 



- + — + 

r 



?)V dhip 



d 



(75) 



r dt J \ dfi 

where ". . ." represents d/d9, d/d(f), and d/dip terms whose explicit form we will not need 
in this spherically symmetric case. 

Finally we consider volume and "surface" elements in phase space. From Eqs. (PH) and 
(^), the spacetime volume element is 

sin 9 



T) 



47rp 



dt AdmAdO A. 



(76) 



Contraction of t] with the coordinate basis vectors d/dt, 8/ dm, 8/ 89, and 8/8(j) produces 
(see Eq. (^)) the spacetime surface elements 

sin 6* 



sin 9 



dm Ad9 A. 



dtAd9 A 



Airp 
sin 6' 

ag = dt A dm A 



Airp 
sin 9 
47rp 
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dt A dm A d9. 



(77) 
(78) 
(79) 
(80) 



We form the volume element on phase space, using Eqs. (|3^) , ([Tq) , and 



Turning to momentum space, from Eqs. (pTD and (|70D-([73D the invariant volume element in 
momentum space is found to be 

Tim = € de A dfi A dip. (81) 

(82) 
(83) 

With the help of Eqs. (|75D and (|83[) , we also construct the surface element in phase space 
that is normal to the phase flow, the 6-form u of Eq. 



sin 6 e 
47rp 



dt A dm A d9 A dcp A de A dfx A dip. 



sin 9 e 



(84) 



[{Lmf dm A d9 A d(j) A de A dp, A dip - (Lm)"^ dt A d9 A d(f) A de A dp, A dip+ 



[LmY dt A d9 A d(f) A dm A dp. A dip - {LmY dt A dm A d9 A d(j) A de A dip + . . .] (85) 

^^^^ ^ dm A d9 A d(f) A de A dp, A dip - sin 9 e^pL dt A d9 A d(j) A de A dp, A dip + 
Anp 

sin 9 



47rp 
sin^ e^(l — /i^ 
47rp 



a / 3f (9 In p 



1 / 3i; 91np 



dt A d9 A dcp A dm A dp, A dip — 

dt A dm A d9 A d(j) A de A dip + 



(86) 



where ". . . " represents terms arising from the d/d9, d/d(f), and d/dip terms of the Liouville 
vector. The validity of Eq. (|37D can be verified (see Eq. (pSf ) for the definition of the exterior 
derivative): 



du 



d ( sin 9 d 
d f sin 6' 



de y ATxp 
d \ sm9 e^(l — p^ 
dp. 



2 , 3t; 91np 



+ 



47rp 



1 f?)V dlnp 

- + — + 

r V r at 



+ 



dt A dm A d9 A d(j) A de A dp, A dip 
0, 



(87) 
(88) 



where the expression for dr/dm is taken from Eq. (|45|). The terms explicitly displayed in 
Eq. (|8^) do in fact sum to zero by themselves, and it can be shown that those represented 
by ". . . " do as well. 
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III. THE DISTRIBUTION FUNCTION, BOLTZMANN EQUATION, AND BAL- 
ANCE EQUATIONS 



In this section — again following Ehlers P2[ for the general case — we define the distribu- 



tion function /, derive the Boltzmann equation, and present the equations describing the 
balance of particle number and 4-momentum. In the general case we reproduce Ehlers' 
derivation of the momentum-integrated number balance equation, and explain why it does 
not yield a conservative reformulation of the Boltzmann equation. We also present the 
balance equation for 4-momentum without derivation. In the special case of Lagrangian 
coordinates in spherical symmetry to 0{v), we verify part of the general derivation of the 
Boltzmann equation by explicit calculation. We also display the Boltzmann equation and 
the momentum-integrated number and 4-momentum balance equations pertaining to this 
special case, including the transformation to the lab frame needed to obtain a "conserved" 
energy. 

A. The distribution function, Boltzmann equation, and balance equations: the 
general case 

The distribution function f for a particle of a given type represents the density of parti- 
cles in phase space. The particle type's 7-dimensional phase space for particles of definite 
mass. Mm, is filled with trajectories (x(A),p(A)), or "states". As a collection of particles 
evolves, the number of particles in each state changes due to collisions. If one considers 
a 6-dimensional hypersurface S in M^, the ensemble-averaged number A^[S] of occupied 
states crossing E is 

iV[S] = Jju, (89) 



where the surface element 00 is given by Eq. (^4]). (Hypersurfaces like S are oriented, and 
the particle trajectories have a direction associated with them as well. The crossing of a 
hypersurface by an occupied trajectory can give a positive or negative contribution to Eq. 
(^), depending on their relative orientation.) With this definition, the distribution function 
f{x,p) is a scalar p2| . 

An equation governing the evolution of the distribution function is obtained by consider- 
ing a closed 6-dimensional hypersurface dD bounding a region D in Mm- The net number 
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of occupied states emerging from D is, from Eq. (^) and the generalized Stokes' theorem, 



N[dD] = fuj= d{fu). (90) 

JdD J D 

(We note in passing — and discuss in more detail in Sec. — that this expression, which 
relates a volume integral to a surface integral, is key to obtaining conservative formulations 
of kinetic theory.) Using the "product rule" of Eq. (P7|), the vanishing exterior derivative of 
uj (Eq. (|37|) ), and the definition of u as the contraction of the Liouville vector Lm with the 
volume element VL (Eq. (H)), we have 

d{fuo) = df ^uJ = df ^{L^■n). (91) 

Any scalar function /, vector field L, and n-form field VL on an n-dimensional manifold obey 
the identity 

df A{L-n) = L[f]n. (92) 
Hence the net number of occupied states emerging from D is 

N[dD] = [ L^[f]n. (93) 
Jd 



(Equations (0)-(93) amount to the Liouville theorem in a relativistic context.) The domain 
D in phase space consists of a region H in spacetime, together with regions Kx in the 
momentum space (the mass shell in the tangent space) at each spacetime point x. Recall 
also that the volume element Q in phase space is the product of the spacetime and momentum 
space volume elements, Q = r] Anm- Hence the integral over D in Eq. (p3|) can be expressed 
as the iterated integral 



N[dD] = v{ L^4f]^m]. (94) 

Because dD is a closed surface, the net number of particles in trajectories emerging from 
D must equal the net number of collisions in D. If correlations between particles can be 
neglected, the spacetime density of collisions can be expressed in terms of a collision integral 
C[/] that depends only on one-particle distribution functions. Therefore, if C[/]7rm 
denotes the spacetime density of collisions, we have 

N[dD]= [ J [ C[/]7r„y (95) 
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Because the regions H and are arbitrary, comparison of Eqs. (0) and (^) shows that 
the evolution of / is determined by 



Lm[f] = C[/], (96) 

or, using Eq. (|2^) , 

p'^'aIt + {eFhp' - rh^p'pn^^ = C[f]. (97) 

This is the Boltzmann equation. 

Next we consider the particle number 4-current, electromagnetic 4-current, and stress- 
energy tensor, and present the balance equations they satisfy. In Eq. (|gD|), specialize the 
hypersurface of Mm to be E = G x P„(x), where G is an infinitesimal spacelike hypersurface 
in spacetime at point x, and Pm{x) is the entire momentum space at point x [^. As 
explained in the endnote to the discussion following Eq. (^6|), on such a hypersurface the 
only relevant term of u is its first term in Eq. (|3^) (given also by Eq. (|36D ), so that Eq. 
(^) becomes 

N[G X Pm{x)] = [cjj [ //:'^^/7r„y (98) 

JG \^Pm(a;) / 

This is the number of particles whose worldlines cross G (at x); hence the 4- vector 

N'-ix) = [ fC-^p^TTm = [ fp'T^m (99) 

JPm{x) JPmix) 

is the particle 4-current density. The electromagnetic 4-current is obtained by multiplying 
by the electric charge: = eN^. The particle 4-current satisfies the balance equation 

and similarly for the electromagnetic 4-current. (Of course the collisions will be such that the 
sum of the divergences of the electromagnetic 4-currents of all particle species will vanish.) 
The stress-energy tensor is 

T^" = [ fC^^C'-.pf'p'nm = [ fp'p'^m, (101) 

which obeys the balance equation 
1 d 



(v^T^^) = F\r - r%,T^^ + J <C[f]C%p%m. (102) 



-g dxf^ 
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To conclude this subsection, we reproduce Ehlers' derivation of Eq. ( p.00|) for number 
balance, and explain why this proof does not yield a conservative reformulation of the 
Boltzmann equation. The proof involves forming an integral over a suitable hypersurface in 
phase space, evaluating that integral in two different ways, and comparing the results. 

First we specialize the integral in Eq. ( P^D to a specific hypersurface of integration. 
Consider an arbitrary region in spacetime, D, whose boundary is dD. Form a hypersurface 
dD in the seven-dimensional phase space M^, composed of the boundary region dlD in 
spacetime together with the entire momentum space Pm{x) at each point of dD. Equation 
(^) becomes 



N[dD] = / fuj, (103) 

JdD 

the integral we will evaluate in two different ways. 

On the one hand, as derived in Eqs. (]90|)-(]93|), the integral in Eq. ( |103| ) can be written 

N[dD] = [ L^[f]n. (104) 
Jd 

With the Boltzmann equation, Eq. (|9 



N[dD] = / Cim = V( C[/]7r„ ) . (105) 

'd Jd \JPm{x) 



On the other hand, the integral in Eq. ( |103|) can be evaluated directly. As a first step, 
integrate over momentum space at each point of dD. As discussed in Subsec. [11 B| , in such 
integrals over 3-dimensional regions in momentum space only the first term of u in Eq. (3l) 



contributes (see in particular the endnote in the sentence following Eq. (^)). Using Eq. 
(PB|) to express the first term of u, Eq. ( |103| ) can be expressed as 



N[dD] = a A jC^f.P^'frfm] (106) 

IdD \Jp,^{x) / 

cr^N'^ix), (107) 

dD 

where the definition (Eq. (|99|) ) of the particle number current has been employed in the 
second step. Equation (|107|) is a closed surface integral in spacetime; by the divergence the- 
orem (a special case of the generalized Stokes theorem), it can be converted into a spacetime 
volume integral: 

N[dD] = [vN",^ (108) 
Jd 

V 4^ii:;;iV^9Nn. (109) 



D V-9^x^' 
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Because the spacetime region ID is arbitrary, Eq. (|100|) follows from comparison of Eqs. 
( |105|) and ([109|). 



Now we can see why this type of proof does not yield a conservative reformulation of 
the Boltzmann equation. In particular, we can see why no direct insight is gained into the 
fate of the non- conservative momentum derivative terms in the Boltzmann equation (Eq. 
(p7|)) upon integration over momentum space. The left-hand side of Eq. ( pTj) arises from 
the action of the Liouville vector on /. Similarly, the phase space surface element uj of Eq. 
(p^) is given by the contraction of the Liouville vector with the phase space volume element 
VL. Inspection shows that the second term of Eq. (0) is closely related to the momentum 
derivative terms in the Boltzmann equation. But in the above derivation of the number 
balance equation, an expression for A^^ is obtained at an early stage by integration over 
momentum space; this integration kills these terms of cu, so that they are no longer present 
by the time the spacetime divergence theorem is applied. The power and elegance of exterior 
calculus has hidden the algebraic details of the connection between the Boltzmann equation 
and the number balance equation! 

This is remedied in Section |V[ Instead of choosing a specialized hypersurface involving 
all of momentum space and evaluating the surface integral of the left-hand side of Eq. (pO|), 
the path to conservative reformulations of the Boltzmann equations involves analysis of the 



volume integral on the right-hand side of Eq. (|90D; for here all the terms of uj have been 
brought inside the exterior derivative, including the ones related to the momentum derivative 
terms appearing in the Boltzmann equation. 

B. The distribution function, Boltzmann equation, and balance equations: La- 
grangian coordinates in spherical symmetry to 0{v) 

The first result to specialize from the previous subsection is the Boltzmann equation. An 
aspect of the derivation that can be demonstrated by explicit calculation is the assertion 

df^uJ = L^[f]n (110) 

contained in Eqs. (PT|)-(P^. In spherical symmetry / = f(t,m,e,fi), so that its gradient is 

df = ^dt + ^;^dm + ^de + ^djj,. (Ill) 
ot dm oe on 
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Forming the wedge product with uj of Eq. (|86D results in 



df Acu 



sin 6 df 



+ r sin 9 e fi — h 



df sin 6 e" 



dm 



Anp dt 
sin 61 e^(l — /i^ 
47rp 

(it A dm A d6 A d(f) A de A dfi A dip, 



ATxp 



3i; 9 In p 
r (9t 



de 



+ 



1 /3f 91np 
r \ r ai 



dp 



X 



(112) 



where all other terms vanish because the the wedge product of a 1-form with itself vanishes 
due to antisymmetry. Comparison with Eqs. ([75|) and (|83|) for Lm and VL then shows that 
df Auj = Lm[f]^, as was to be demonstrated. The Boltzmann equation, Eq. (|97D , specializes 
to 

df 



df . 2 df 2 



2 I d\np 



+ 



6(1 -p^ 



1 /3f (91np 

- + — + 

r V r at 



dp 



C[/]. 



(113) 



This agrees with Eq. (20) of Ref. |2| . 

Next we specialize Eq. ( |100| ) for particle number balance. Equations (0), (|5lD-(p4D, 
(^)-(|7^), and (|STD, together with the fact that / = f(t,m,e,p) in spherical symmetry, 
imply that 



AT* 



f 271 de dp = pH 



N 



Airr^p I f 2-ne^pdp = A'nr'^p^G^ 



(114) 
(115) 



are the only nonvanishing components of the particle number vector. With these expressions 
and Eqs. <^ and (^, Eq. (^) becomes 

dH^ 

dt am ^ ' ' p 



+ — (ATxr'^pG^) = - / C[/] ededpd^. 
dm ^ p J 



(116) 



This agrees with Eq. (28) (with lapse function a set to 1 for flat spacetime) of Ref. [|T9[] . 

Finally we specialize Eq. (|102|) for particle number 4-momentum balance. Equations 
(H), (|51|)-(|D, (^-(^), and (|lD, together with the fact that / = /(t, m, e, p) in spherical 
symmetry, imply that 



T" = J f 2Tre^dedp = pH^, 
jtm^jrut ^ 4^^2^ f f 2Tre^pdp = A7rr^p^G^, 



(117) 
(118) 
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j^mm ^ (^4T^r^pf J f 2TTe^^'^d^= {Anr^pfpF^, (119) 
T'' = y,j f vre^ll - /i^) dp ^£r2iH^- P"")^ (120) 

T'^^ = / / ^^'(1 - ^^') df^ = TTY^iH^ - P"") (121) 

sm 6 J Ir^ sm d 

are the only nonvanishing components of the particle stress-energy tensor. With these 
expressions and Eqs. (^8]) , (l55| ) - (l65|) , (0)-(0), (|70D-([7^), and (|8T|) , the t and m components 
of Eq. ([T02|) become 

am ^ ^ \r at ) 

C[/] dedpd^Xlll) 



dt 



1 9 



With the help of Eq. (|45|), the m component equation can be expressed 



dt 



= - / C[/] e^pdedpdifi. (124) 

Equations (|122|) and (|124|) agree with the expressions following Eq. (31) of Ref. [|l^. As 
expected from the discussion in Sec. |, these 4-momentum balance equations have non- 
conservative terms arising from the non- vanishing connection coefficients of Eqs. (^)-(^). 

It was also noted in Sec. || that in Eulerian coordinates in flat spacetime, the connection 
coefficients T^^i, vanish, so that a conserved energy can be defined. Equation ( p.Q2| ) for 
particle 4-momentum balance is a vector equation; the t component of a vector V is obtained 
from the t and m components by the transformation 

where the transformation coefficients are given by Eq. (^). Equation ( |122| ) is actually p~^ 
times the t component of 4-momentum balance, and Eq. (3^) is [Airr'^p'^)^^ times the m 
component. These factors must be restored before plugging into Eq. ( |125| ); the final result 
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for the t component of 4-momentum balance is 



- I C[f] {l + vfi)€^dedfid^. (126) 



(The 0{v) baryon conservation expression 



^ + — + A-nr^p— = 0, 127) 
ot r om 

obtained from Eqs. (^Op and (^)-(^), has been employed.) Equation ( [126| ) agrees with 
Eq. (32) of Ref. |[T^. The non-conservative third term on the left-hand side is O(f^) and 
ought to be dropped at the level of approximation we are using. However, it is retained in 
Ref. |jT9| for a practical reason. In supernova simulations velocities can exceed values for 
which the 0{v) approximation is strictly valid; and while use of the 0{v) formalism might 
be questioned on physical grounds, the O(w^) non-conservative term can still be used to 
check that the numerical implementation of the 0{v) Boltzmann equation is consistent with 
Eq. (11261) . 



IV. CONSERVATIVE FORMULATIONS OF PARTICLE KINETICS 



We are now in a position to present conservative formulations of kinetic theory. We seek 
expressions closely tied to Eq. (|100|) for particle number balance and Eq. (|102|) for particle 
4-momentum balance. In addition to deriving these conservative formulations in the general 
case, we show the relationship between them and explain the care that must be exercised in 
finite differencing the number balance equation in order to make it consistent with the energy 
balance equation. We also specialize these results to Lagrangian coordinates in spherical 
symmetry to 0(f), and see that the "number conservative" Boltzmann equation — arrived 
at in the past by guesswork — emerges naturally from our formalism. 



A. Conservative formulations of particle kinetics: the general case 

For an expression related to particle number balance, the derivation of the Boltzmann 
equation in the previous section had us closer to the desired result than might first be 
realized. In rushing headlong towards an equation for /, the key relation is easily overlooked: 
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It is Eq. (0), the result of the generahzed Stokes' theorem. The integrand on the right-hand 
side of Eq. (pOD, d{fuj), is conservative: Having been obtained from Stokes' theorem (the 
generahzation) of the divergence theorem, it has everything inside the exterior derivative. 
Being an exterior derivative, it is too abstract to be directly useful; but massaging it just 
enough to bring it into the form d{fuj) = N[f]Q, where Q is the volume element in phase 
space, we shall see that N[/] is in fact a conservative differential operator in the familiar, 
elementary sense. 

Now we take a detailed look at Eq. (|90D , in particular the exterior derivative 



difu;) 



d{^jfC^f,fe,\,,„\dx''n A vr™ - 
/ 1 

d 



E(p) 
1 



det^ 
au 



du''^ ) A d7], 



;i28) 



where we have used Eqs. (0), (pS]), (p6D, and ( |27D and employed the vertical bar notation 
introduced in Eq. (p^). The exterior derivatives will be expressed in terms of the basis 
{dx^,du^} on Mm- First we note that in our chosen momentum coordinates, the second 
term in Eq. ( 128 ) vanishes because diim = 0: This is because Hm as expressed in Eq. (^Tj) 
has no dependence on {x^}; and while it depends on adding another momentum 1- 

form to the wedge product du^'^^ would cause it to vanish due to its antisymmetry. Similarly, 
the fourth term vanishes because drj = 0. To understand what happens to the first term, 
consider one of the terms in the sum over fi: 

d , 



123\ 



A TT^ 



gp''C\f)dx 



0123 



A 7r„ 



(129) 



While derivatives with respect to {x*} and are nonvanishing, only the wedge product 
of dx^ with dx^"^^ is non-vanishing, and the wedge product in does not admit another 
momentum 1-form. The other terms in the sum over /i are similar. To understand the third 



term in Eq. (p.28|), consider a particular term in the sum over i: 



d 



1 



E(p) 

d 



det 



dp. 



du^^ ) A r] 



dph> 



du^ \E{p) 



det 



[-] 



01 23 



icF^op" - T^.^p-^pP)-^/ ]du''' A V 
op^ 



(130) 
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and similarly for the other terms in the sum over i. With Eq. (|33|) , all terms in Eq. (|128|) 
can be assembled to form the expression 



difcu) = N[f]n, 



:i311 



where 



1 



d 



det|§P 
ou 



1 



d 



1 



op) 



(132) 



In the derivation of the Boltzmann equation, Eq. (|131|) can be used to replace Eq. (|93| 
with a similar expression in which L^if] is replaced by N[/], with the result that 



N[/] = C[/]. 



(133) 



This is our conservative formulation of particle number kinetics. 

Next we seek a formulation related to particle 4-momentum balance. With arbitrary u^, 
we evaluate the exterior derivative d{v^p^fuj) in two separate ways. (To simplify formulae, 
we often write p^ instead of C'^jxp'^, but because of our chosen momentum coordinates on 
Mm the latter expression must often be used in computational steps.) First, employing Eqs. 
(HI) and (|131), 

d{v^p^fuj) = d{v^p^) Afuj + v^p^C[f]n. (134) 

Computation shows that 



d{v^pn A fu 



,1 



(135) 



Because of the mass shell constraint (Eq. ([T9|)), p° is considered a function of the {p*}. The 
geodesic equation (equation (|1^)) can be used to show that 



^{eFhp' - Thpp'p') = eFKp' - T^pp'^pK 
op^ 



(136) 

This, together with Eq. (|T^) for the transformation of the connection coefficients, allows 
Eq. (|135|) to be written 



d{v,pn Afuj= {p'p'f^ + v,{eF^up' - r^,pV)/)fi. 



(137) 
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Second, we evaluate d{v^p^fuj) directly. A computation similar to that leading to Eqs. 
(^)-(|T32l) yields 



;i38) 



where 



1 d 



-g dx^ 



-gC%C%pl^p' f) + 



au 



-1 



d ( 1 



au 



{eF\p^ - T^,-p^pP)-^C>^-p^f ) .(139) 



Recalling that is arbitrary, Eqs. (|134|) and (|137|) - (|139|) can be combined into 



(140) 



a "conservative" formulation of particle 4-momentum kinetics. 

To illuminate the practical problem of accurately accounting for both particle number 
and energy, we consider the relationship between Eq. (|133| ) for particle number kinetics and 
the t component of Eq. ( |140D for particle 4-momentum kinetics. (This is an issue when 
one follows the general approach of, for example, Refs. ^, Q , in which a particle number 
distribution equation is solved and a particle energy equation is used as a consistency check.) 
Specifically, we need to relate the spatial and momentum divergence terms of these equations 
to each other and identify terms that cancel. 

First we relate the spatial and momentum divergence terms of Eq. ( |133|) and the t 
component of Eq. (|140|) . Their spatial divergence terms, 

1 d 



1 

'—g dx^ 



are related by 



-gC^P^f), 



^gC^C^i^p^p^f) 



d 



Their momentum divergence terms 
(T*)^[/] = E(p 



det 
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[-] 


det 


[^] 
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(141) 
(142) 

(143) 
(144) 



- - ~ (9?/* 

{eF^-p>^-V^^p^pP) f 

op> 

- - (9?/^ - \ 
(eF^^p^ - P^^pV)^/:*PPVj;45) 



29 



are related by 



(146) 



Comparison of Eqs. ( |143| ) and ( |146| ) with Eq. ( |133| ) and the t component of Eq. ( |140| ) 
shows that the following equation must be valid: 



Es + Em = eF'^C'^f.pf^ - T\pC^ (,C%p^p^ 



(147) 



where 



:i48) 



is the "extra" term of Eq. (|143|) relating the spatial divergence terms of the number and 
energy balance equations, and 



Em = f {eF^f^P^ - T^'^ppV)^ ^^'^ 



(149) 



is the "extra" term in Eq. (|146|) relating the momentum divergence terms of the number 
and energy balance equations. Noting that 



^ ^ = + s'^c'- 
dpi " du" dp^ ° i 

Eq. ( |136| ) can be used to rewrite Eq. ( |149| ) as 

Use of Eq. (|T7D and two applications of the identity 
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lead finally to 
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Equations ( |148| ) and (|153|) make it obvious that Eq. (|147|) is satisfied. 

In a computational approach to transport in which a particle number distribution equa- 
tion is solved and a particle energy equation is used as a consistency check, care must be taken 
that Eq. (|147|) is satisfied numerically. In particular, a given finite difference representation 
of the spacetime divergence ^x[f] implies a corresponding finite difference representation 
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of position-dependent quantities in Es, and the position-dependent quantities in terms of 
Em that cancel with Es must have a matching finite difference representation. Similarly, 
a given finite difference representation of the momentum divergence Np[/] implies a corre- 
sponding finite difference representation of momentum variables in Em, and the momentum 
variables in terms of Es that cancel with terms of Em must have a matching finite difference 
representation. 



B. Conservative formulations of particle kinetics: Lagrangian coordinates in 
spherical symmetry to 0{v) 



Using Eqs. (||), (|T])-0, ( 
tion of Eq. ( |133| ) specializes to 

p J am \ p 



and (|7DD-(|73D, the "number conservative" formula- 
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which agrees with Eq. (23) of Ref. ||2^. In Ref. the necessity of making f/p the 

evolved variable in order to get a conservative particle number equation is left unexplained. 
In the present formalism, we see that the factor 1/p comes from the factor y/—g = sinO/Airp 
in Eq. (|100|) . Multiplication of Eq. ( |154| ) by dedpdip and integrating immediately yields 
Eq. (|116|) , expressing particle number "conservation". 

Similarly, the t and m components of the "momentum conservative" formulation of Eq. 
[40|) are 
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and 
1 d f r'^pepf 



r'^p dt 
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^ f2v d\np\ ^ 2 , / 1 
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(156) 



With the help of Eq. 



the m component equation can be expressed 
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Multiphcation of Eqs. (|155|) -( [l57| ) by e^dedpdip and integrating immediately yields Eqs. 
( |122|) - (|124|) , expressions for particle 4- momentum balance in the comoving frame. 

Just as the comoving frame t and m momentum-integrated balance equations — Eqs. (|122|) 
and ( |124|) — can be combined to form the lab frame t "conservation" equation, Eq. (|126D , 
Eqs. ( |155| ) and (|157|) can be combined as 
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Multiplication of this equation by de dp dip and integrating immediately yields Eq. (|126|) 
for lab frame energy "conservation" . 

We now examine the relationship between Eq. ( |154| ), a conservative formulation of par- 
ticle number kinetics, and Eq. ( |158| ), a conservative formulation of particle energy kinetics. 
We multiply Eq. ( |154{ ) by (1 + vp) and consider what it takes to get the terms in Eq. ([1581 ). 
In doing so, we examine carefully only 0{v) terms. First the time derivative terms: 
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Next we compare the mass derivative terms: 
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In the notation of the previous subsection, the "extra terms" from the spacetime divergence 
are 

. r)i) f 

(161) 



E5 = 47rrV|^- + 0(t;2). 



dm p 

Next we relate the momentum divergence terms in the number and energy equations, looking 
first at the energy derivatives. 
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and then at the angle derivatives: 
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P/ r p 

In the notation of the previous subsection, the "extra terms" from the momentum divergence 
are 
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Having Es and Em, we are ready to verify Eq. ( [L47| ), which must be satisfied for consistency 
between the particle number and energy balance equations. Having specified electrically 
neutral particles, and having chosen to work with the lab frame energy expression of Eq. 
( |158|) (in which the connection coefficients vanish), the right-hand side of Eq. ( |147| ) vanishes. 
Employing Eq. ( |127| ) for baryon conservation, from Eqs. ( |161|) and (|164|) it is easy to see 
that 

Es + Em = (165) 

is indeed satisfied analytically. But in solving the "number conservative" Eq. (|154|) numeri- 
cally, consistency with with the "energy conservative" Eq. (|158|) requires that Eq. (|165|) be 
satisfied numerically as well. Ref. [H provides an example of a finite difference representation 



of Eq. ( 154 ) that satisfies this criterion. 
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V. CONCLUSION 



In this section we summarize our conservative formulations of kinetic theory, comment 
on their relation to moment formalisms, and discuss their possible application in the core- 
collapse supernova environment. 

Having in mind computational radiative transfer in astrophysical environments, we have 
sought formulations of relativistic kinetic theory with the following properties: (1) they 
are expressed in terms of global, Eulerian (or "lab-frame") spacetime coordinates {x^}; (2) 
they are expressed in terms of convenient three-momentum coordinates (e.g. spherical 
polar), which are taken from the orthonormal momentum components {p*} measured by an 
observer comoving with the medium; and (3) they are "conservative", having transparent 
connections to total particle number and 4-momentum balance as expressed in Eqs. (|100|) 



and ([101 . 

To express our formulations having these properties, we here introduce the specific particle 
number flux vector 

= C'^f.pf'f (166) 
and the specific particle stress-energy tensor 



(167) 



where the transformation to the comoving frame ^C^^ is given by equation (|I^). (While the 
adjective "specific" often denotes a quantity measured per unit mass, in this context we 
use it to denote the particle flux and stress-energy in a given invariant momentum space 
volume element.) While the distribution function / of a given particle type of mass m and 



charge e obeys the Boltzmann equation (Eq. (|97|) ), the specific particle number flux and 
stress-energy satisfy the conservative equations 
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p^^^^u _ Y^^^^r^p + /:^^/c[/] . (169) 



In stating that Eqs. ( |168|) and (|169|) constitute conservative formulations of particle kinetics, 



we mean that the connection to the balance equations of Eqs. ( |100D and (|102|) is transpar- 
ent, in the following sense. We can use elementary calculus to form the familiar invariant 
momentum space volume element 



d^p 1 



E{p) E{p) 



ou 



d\ (170) 



where a transformation from orthonormal momentum components {p*} to some other set 
of coordinates (e.g. momentum space spherical coordinates {u*} = {|p|, -i?, </?}) has been 
performed. Multiplying Eqs. ( |168| ) and ( |169| ) by Eq. ( |170| ) and integrating, the terms with 
momentum space derivatives are obviously transformed into vanishing surface terms; the 
results are Eqs. ( |100| ) and ( |102| ) for total particle number and 4-momentum balance. 

In terms of differential forms, the procedure for obtaining the conservative formulations 
of kinetic theory is straightforward. First, express the volume element Q in the phase space 
for particles of definite mass in terms of the desired spacetime and 3-momentum coordinates. 
Next, by contraction with the Liouville vector, form the hypersurface element uj = ■ Q. 
Then bring the exterior derivative d{fuj) into the form by direct computation; Eq. 

( |168|) results on comparison with the Boltzmann equation. This result can then be used in 
conjunction with an evaluation of d{v^p^fuj) for arbitrary to obtain Eq. ( [L69| ). The reason 



the procedure is straightforward is that the "heavy lifting" of transforming the Boltzmann 
equation into conservative forms is handled by two "levers" of considerable power: the 
generalized Stokes theorem, and the key relation rfa; = 0, which is closely related to the 
relativistic Liouville theorem. 

A concrete example of our formalism is provided in the appendix, which contains 0{v) 
equations for the specific particle number density, specific particle energy density, and their 
angular moments — all in flat spacetime, but in coordinates sufficiently general to represent 
rectangular, spherical, and cylindrical coordinate systems. 

We now comment on the connection of Eqs. ( |168| ) and (|169|) to moment formalisms. In 
the usual treatments, if one writes the distribution function as a function of momentum 
variables as measured in a given frame (lab or comoving), it is natural to form moments 
by multiplying the distribution function by, for example, energies and angles measured in 
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that frame, and integrating. Lo and behold, it turns out that these moments are number 
densities and fluxes, and energy /momentum densities and fluxes: components of a particle 
number flux vector and stress-energy tensor, measured in the same frame chosen to measure 
angles and energies. Traditional, then, are treatments in which the components of conserved 
tensors as measured in a given frame are expressed as functions of momentum variables as 
measured in that same frame. 

But this traditional approach to moments may not be the most convenient, and Eqs. (|168|) 
and ( |169| ) provide an attractive alternative. For example, Liebendorfer et al. |19| form 
moments in the traditional way, resulting in components of conserved tensors as measured 
in an orthonormal comoving frame. But it is the tensor components in the lab frame that 
one would like to check, either because it is the coordinate basis (at least in a spatially 
multidimensional simulation) and therefore natural to deal with, or because it is the basis in 
which energy is conserved (in a simulation in comoving coordinates in spherical symmetry). 
The required transformation of the tensor components between these frames leads to the 
numerical complexity mentioned in Sec. |, and discussed further in the latter parts of 
Subsecs. IV A and [IV B| . In contrast, when the coordinate basis is in the lab frame as 
expected in spatially multidimensional simulations, integration of Eqs. I \168i ) and l \169{ ) over 
momentum variables leads directly to the tensor components in the desired coordinate basis, 
even though the specific particle number and stress- energy are functions of comoving frame 
momentum variables. The insight here is that the frame in which the tensor components are 
measured need not be the same as that employed to obtain the momentum variables used 
to parametrize the particle distributions. The appendix provides an example of a moment 
formalism of this kind. 

In simulations of systems like core-collapse supernovae — in which careful attention to en- 
ergetics is critical — a number of possible approaches, based on the conservative formulations 
of kinetic theory presented in this paper, might be suggested. 

First, the general approach of, for example, Refs. p, |^, |[| could be followed, in which the 
conservative particle number distribution equation (Eq. (|168|) ) is solved, and the conserva- 
tive particle energy equation (the time component of Eq. ( |169| )) is used as a consistency 
check. The quantitiy y/—gAf^, which might be called the specific particle number density, 
would be the primary neutrino distribution variable: It is the contribution of each comoving 
frame momentum bin to the lab frame particle number density. This approach makes num- 



36 



ber conservation a somewhat natural outcome, but energy conservation would require finite- 
differenced representations of various quantities to be "matched" in order that Eq. ( |147| ) 
be satisified numerically. This might be considered the most rigorous and self-consistent 
method. (Note that if one solves the "plain", non-conservative Boltzmann equation — Eq. 
(P7|) — for the scalar distribution function / as a function of comoving frame momentum vari- 
ables, conservation of neither lab frame particle number or energy is straightforward. The 
same is true of methods (e.g., Ref. [^) based on a non- conservative form of the transport 
equation for the comoving frame specific intensity.) 

In order to avoid the intricate finite differencing of numerous terms demanded by this 
method, a second option would be the use of T"", which might be called the specific par- 
ticle energy density, as the primary neutrino distribution variable. Designing a code around 
a differenced version of the fi = component of Eq. (|169|) would make accurate accounting 
of total neutrino energy (as represented by the /i = component of Eq. ( |102| )) relatively 
straightfoward. Of course, the neutrino number balance equation expressed in terms of 
y/—g would have numerical errors unless certain finite differencings were carefully de- 
signed. But with respect to the crucial energetics of the physical system, it is worth noting 
that there are a couple of factors mitigating the impact of errors in number conservation in 
comparison with errors in energy conservation. Errors in number conservation translate into 
errors in the electron fraction Yg, which affect energy conservation through the equation of 
state, but: (1) Only and (z/)^ affect Y^, while all species impact the energy budget. Better 
to have two species contributing to error rather than six! (2) The effects of z/g and (u)^ on Ye 
are opposite in sign (unlike their contributions to energy), so that to the extent that their 
distributions are similar, the impact of their errors on Y^ may approximately cancel. 

A third possibility would be to solve for both the specific particle energy density y/—g 
and the specific particle number density ^J—gM^. (Rampp and Janka [^] solve for both 
number and energy distributions, but in the comoving frame; this limits the utility of their 
approach with respect to accurate tracking of lab frame quantities.) Instead of pre-defining 
both the boundaries and center values of bins in energy space, one could define the bound- 
aries only and use the values of \/—g T^^ and \/—gM^ (along with the transformation C^^) 
to obtain center values of the energy bins in each spatial zone and each time step. The con- 
sistency of the solutions would be arguably reasonable as long as the derived center values 
of the energy bins do not wander outside the pre-defined bin boundaries. 
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APPENDIX: FLAT SPACETIME, 0{v) EQUATIONS IN A GENERAL COORDI- 
NATE SYSTEM 



In this appendix we present an application of the conservative formulations of kinetic 
theory derived in this paper: 0{v) equations for the specific particle number density and 
specific particle energy density, as well as angular integrals of these equations — all in flat 
spacetime, but in coordinates sufficiently general to represent rectangular, spherical, and 
cylindrical coordinate systems. The angle-integrated equations constitute "monochromatic" 
moment formalisms that provide an alternative to traditional variable Eddington factor 
methods of handling radiation transport. We specialize to massless, electrically neutral 
particles. 

We begin by describing our spacetime and momentum space coordinate systems. While 
we assume fiat spacetime, in order to accomodate curvilinear coordinate systems we employ 
a general spacetime coordinate labeling (x^) = {x^,x'^,x^,t)^ . The line element is 



ds^ = g^^dx^dx'^, 



(A.l) 



with 



i9, 



(A.2) 



/I Q \ 

a?{x^) 

62(^1)^2(^2) 

\0 -1/ 

In this matrix expression for the metric components, rows and columns are ordered 

1, 2, 3, 0. In Cartesian coordinates, (x^, x^, x^) = {x, y, z), and (a, b, c) = (1, 1, 1). In spheri- 



cal coordinates, ( 



,6,(/)), and {a,b,c) = (r, r, sin^). In cylindrical coordinates. 



{x^,x^, x^) 



{r,z,(l)), and (a, 6, c) = (l,r, 1). Our spacetime coordinate systems are "lab 
frames," so that the equations we derive are Eulerian. Orthonormal "lab frame" coordi- 
nates, indicated by barred indices, are obtained by the transformation 



with 



dxf" = ef'^dx", 

/I 0\ 
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VO 1/ 
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(A.4) 



38 



The transformation to an orthonormal frame comoving with the fluid, indicated by indices 
adorned with a hat Q, is 

dx^ = A'^-^dx'', (A.5) 

where to 0{v) 



(A^ 
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-vi\ 
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-V2 




1 / 



(A.6) 



The bars on the velocity indices are a reminder that the fluid velocity is expressed in the 
orthonormal lab frame coordinate system. The combined transformation from the lab coor- 
dinate frame to the orthonormal comoving frame is 



where 



(A.7) 



(A.8) 



The neutrino 4-momentum is described in terms of its comoving frame components, 
pi- ' 



p^,p'^,p^,p^ ) . Only three momentum variables are independent; we choose polar 



coordinates in momentum space, defined by 



p^ = ecos-i?, 

p^ = e sin'i? cos(/9. 



(A.9) 
(A.IO) 
(A.ll) 



The radiation field is a function of the variables t J Ob J Ob J Ob , e, "i?, (yj. The invariant spacetime 
volume element in Eq. (|28|) is ahcdx^dx^dx^dt, and the invariant volume element on the 
mass shell in momentum space in Eq. (^) is esm-d dedd dip. 

We first present an equation for the specific particle density A/", defined by 



M = p'f = C'f^p^f, 



(A.12) 



where / is the invariant particle distribution function defined in Sec. |lTl| . The specific 
particle density is related to the lab frame particle density n = A^* by 



Afe sin.'d de dd dip. 
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(A.13) 



From Eq. ( |16^ ), we find that the specific particle density satisfies 



Bt[M] + D^i [TV] + D^2[7V] + D^3[7V] + ©^[AT] + ©^[AT] + OJ7V] + 0^[7V] + 0^[7V] = C, 

(A.14) 

where C is the invariant collision integral appearing in the Boltzmann equation (Eq. (P7D), 
and 

DJAT] = ^, (A.15) 
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1 /3 22 3 23\ 



^^(-4C¥:.-S^t^i - 4s^t^2 - Ctf44^2 + 0^4^9.4^3))-^) • (A.23) 

In these expressions we have employed the notation = sin.'d, = cos??. In accordance 
with the 0(f) limit, we have ignored acceleration terms and terms nonlinear in velocities 
and/or velocity spatial derivatives. 

Next we consider an equation for the specific particle energy density S, defined by 

S=pYf = C'f,C'opf'p'f. (A.24) 

The specific particle energy density is related to the lab frame particle energy density e = T** 
by 

e = j £ esin^ded'ddifi. (A.25) 
From Eq. ( |169| ), we find that the specific particle density satisfies 

Dt[^]+D^i[^]+D^2[f]+D^3[^]+D^[^]+D^[f]+0,[f]+0^^[^]+0^[f] = eil+vic^+v-2S^c^+v-sS^s 

(A.26) 

The terms on the left-hand side are given by Eqs. ([A.15| )- ([A.23| ) with A/" replaced by S.^ 



The t component of the second term of Eq. ( |169| ), which we refer to as the "acceleration 
term", vanishes for the fiat spacetime metric of Eq. (|A.2| ). 



We now consider angle-integrated versions of Eqs. ( |A.14| ) and ( |A.26| ), which constitute 



equations for "monochromatic moments" of the radiation field. These moments are functions 
of the variables t, , e. In deriving and presenting these expressions, it is convenient 

to define 

nj = cos^9, (A. 27) 

= Sint? COS(y9, (A. 28) 

ng = sin-i^sinyj, (A. 29) 
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the orthonormal comoving frame components of the unit 3- vector specifying the particle di- 
rection. In accordance with the 0{v) hmit, we ignore acceleration terms and terms nonlinear 
in velocities and/or velocity spatial derivatives. 

Om^ first equation of this type dictates the evolution of the lab frame monochromatic 
particle density defined by 

J Msm'd (I'd dip. (A.30) 
It is related to the lab frame particle density n = A^* by 

n = J^ede, (A.31) 

and its evolution is governed by 

Bt [m\ + D^i m + D^2 [m] + D^3 [m] +Oe[% = J C sin^? d'd dip, (A.32) 

where 

= ^, (A.33) 

= + - ^ii^l - - /i3^3)^), (A.34) 

= ^c'i^W^ + V-2- hiV-l - - /23^3)^), (A.35) 

= ^^((/3 + ^3 - fisV-l - /23^2 - hmm, (A.36) 

1 da . „ „ , 

-^(/22-l-/i2-2) + 
5^(/33^^1-/i3^3) + 

In these expressions, we have defined 

fi = J^n-sm^d-ddcp, 

fl] = Hn-n-sm^d^d^, (A.38) 
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which are akin to Eddington factors in traditional moment approaches to radiation transport. 

Finally, we present a set of coupled equations for the lab frame monochromatic particle 
energy density € and orthonormal lab frame monochromatic particle momentum density 
defined by 

€ = J Ssm^d^dip, (A.39) 
qj^ = e\ J fpYsm^d^dcp. (A.40) 

They are related to the lab frame particle energy density e = T*^ and orthonormal lab frame 
particle momentum density = e'^T*^ by 



eede, (A.41) 

^'ede. (A.42) 

The evolution of the lab frame monochromatic particle energy density is governed by 

B>t[e] + D^i[€] + B>^2[(a] + D^3[€] + Oe[(B\ J C{1 + vini + V2n^ + Vsn^)sm'&di9d(p, 

(A.43) 
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-9^(%^2-%^3))^)- (A.48) 
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The evolution of the first component of the orthonormal lab frame monochromatic particle 
momentum density is governed by 



(A.49) 



where 



mV'] = (A.50) 

D.i[^'] = ^-^{ab{2vi^^ + {hii - 2{hrnvi + hi^v^ + h^.m))'^)): (A.51) 

D,2[^I] = lA.(^c{^% + ^^y-^ + ^hii-2{hr,m + k22V-2 + h^ (A.52) 

Dx3[^^] = ^^(^'^i + ^'^3 + ihs - 2(/iii3^^i + h^^V2 + h^^v-s))e), (A.53) 

F22[q3^] = ---^{2^% + - 2(/li22^1 + %2^2 + /i223^3))«), (A.54) 

F33[^'] = ~-g^{2^% + (% - 2(/ii33^;i + %3^2 + %3^3))^), (A.55) 
^ r^Ti 1 <9 / 9 / , 9^2 , (5^3 , 
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]_f^h M^h \m 

bc\dx^ ''^^ dx^ '^^^ dx^ '^y^ 

-^(^122^1-^112^2) + 

1 dc \ \ 

— ^(^i33^'2-/ii23^3)j^j- (A.56) 

The evolution of the second component of the orthonormal lab frame monochromatic particle 
momentum density is governed by 

(A.57) 

where 

nV'] = (A.58) 
= ^^(a6(^^^i + ^S + (/ii2-2(/iii2t'i + /ii22^2 + ^^^^^ (A.59) 
D,2[^2] = 1 A(c(2q32^2 + (hi - 2(/ii22^i + ^2^2 + h^^vsM)), (A.60) 
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I + 
+ 



D,3[^2] ^ i._|_(^3^_ ^ ^2^_ ^ _ 2(/ij23^j + ^3^2 + %3^3))^), (A.61) 

Fl2[^'] = -^(^'^1 + ^'^2 + {hi2 - 2(/lii2^^I + /li22^2 + /ii23^3))«), (A.62) 

1 9c 

F33[r] = --^(2^\3 + (% - 2(/ii33^i + %3^^2 + %3^3))^), (A.63) 
L(^h M^h M^h W 

6cUa;3''i23 + 5^3/^223 + 9^3 f^m) + 

If^h M^h -^^h \^ 

a \dx^ '^^"^ ^ 9x2^222 + 5^2^223; + 

~^^(^222'^1 ~ ^i22'"2) 
"^^(%3^1 ~ ^123^3) 
^^(^233^2 -/i223^3))«)- (A.64) 

The evolution of the third component of the orthonormal lab frame monochromatic particle 
momentum density is governed by 

(A.65) 

where 

= (A.66) 

^A'f] = ^^(«K^'^'i + ^S + (/^i3-2(/iii3t)i + /list's ^ (A.67) 

D,2[^3] ^ + ^'^3 + (% - 2(/li23^1 + him + %3^3))^)), (A.68) 

»x3[^^] ^ ^^(2^3^3 + (% - 2(/ii33^i + %3^2 + %3^3))e), (A.69) 

1 d}) 

Fi3[^'] = -^^(^'^1 + ^'^3 + {h^ - 2(/iii3^i + him + ^133^3))^), (A.70) 
1 dc - 

F23[^'] = -^("^'^2 + ^'t^3 + - 2(/ii23^I + ^3^2 + ^a^s))^), (A.71) 

1/"^^ w 

a \dx^ ^ dx^ ^^^ ^ dx-' ^^V + 

L(^h M^h \m 
be \dx^ ^^3 + ^^3 ^^233 + ^^3 ^^333; + 

-^fe^I- ^123^2) + 
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^^(%3^1 - ^133^3) + 

1 dc \ \ 

-^fe2-/i233^3)jej. (A.72) 

In these energy and momentum equations, we have defined 

1 f 

h-fj = — £ ninjsm{}d{} dip, (A. 73) 

which are akin to Eddington factors in traditional moment approaches to radiation transport. 
We point out that /• is not equal to ^V^5 ^^at the second angular moment factors 
f~fj and /?-•-• are not equal to each other. This is because Af and S have different angular 
dependencies, as can be seen from Eqs. (|A.12| ) and ( |A.24| ). This complication is a result 
of taking moments of lab frame number and energy distributions with respect to comoving 
frame angular variables. 
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In a "coordinate basis" the basis vectors are = ^|jr • These basis vectors "commute" because 
the order of partial derivatives can be interchanged freely. The connection coefficients in a 
coordinate basis are given by V^j^p = \g^'^ {^§§^ + — ^§^^- In a "non-coordinate basis" 
obtained from the coordinate basis by a transformation e^/ = ^/e^, the basis vectors e^i = 
^^m' af^ do not commute if L^^i depends on position in spacetime. The connection coefficients 
in a non-coordinate basis have extra terms associated with the non-vanishing commutators 
of the basis vectors. In this paper, the connection coefficients in a non-coordinate basis are 
obtained from those of a coordinate basis via the transformation T^^' yipi = L^^' ^^L^ yi piT^^ yp + 

The contraction of a basis vector with a 1-form dz^ is ■ dz^ = <^a- In the contraction 
of a vector with a direct product of 1-forms, the contraction is with the first 1-form in the 
direct product. This means that the contraction of a vector with a wedge product of n 1-forms 
gives n terms, because the wedge product is the completely antisymmetrized sum of n direct 
products. 

The computational procedure for integrals over differential forms involves (1) forming an in- 
finitesimal parallelpiped — a wedge product of the vectors forming the edges of the infinitesimal 
region — at each point in the region over which the integration is performed; (2) contracting 
this infinitesimal parallelpiped with the differential form at each point; (3) calculating the in- 
tegral with the usual rules of elementary calculus. An infinitesimal parallelpiped in momentum 
space is a wedge product of displacement vectors in each of the three directions in momentum 
space. In the contraction of this parallelpiped with lo, the only terms that do not vanish are 



those in Eq. (pq); the other terms of Eq. (34) have only two momentum space 1-forms. 
Other 6-dimensional hypersurfaces can be formed, which contain a 4-dimensional region in 
spacetime and a 2-dimensional surface in momentum space or a timelike region in spacetime 
together with a 3-dimensional region in momentum space; such hypersurfaces are not of use 
in defining the particle number 4-current and stress-energy tensor. 

Because the relations between J\f , E and / have different dependencies on velocity, expressions 
of Eqs. (A. 14) and ( |A.26| ) in terms of / would differ by terms of O(u^) and higher. While 



in principle these discrepancies vanish to 0{v\ a simulation solving for both M and E in 
which turns out to be non-negligible compared to unity can be expected to exhibit some 
inconsistency between values of / derived from Eq. (A. 12) and Eq. (A. 24). 
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